A Two-Region Diffusion Model for Current-Induced Instabilities of Step Patterns on 
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We study current-induced step bunching and wandering instabiUties with subsequent pattern 
formations on vicinal surfaces. A novel two-region diffusion model is developed, where we assume 
that there are different diffusion rates on terraces and in a small region around a step, generally 
arising from local differences in surface reconstruction. We determine the steady state solutions for 
a uniform train of straight steps, from which step bunching and in-phase wandering instabilities 
are deduced. The physically suggestive parameters of the two-region model are then mapped to 
the effective parameters in the usual sharp step models. Interestingly, a negative kinetic coefflcient 
results when the diffusion in the step region is faster than on terraces. A consistent physical picture 
of current-induced instabilities on Si(lll) is suggested based on the results of linear stability analysis. 
In this picture the step wandering instability is driven by step edge diffusion and is not of the Mullins- 
Sekerka type. Step bunching and wandering patterns at longer times are determined numerically 
by solving a set of coupled equations relating the velocity of a step to local properties of the step 
and its neighbors. We use a geometric representation of the step to derive a nonlinear evolution 
equation describing step wandering, which can explain experimental results where the peaks of the 
wandering steps align with the direction of the driving field. 
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I. INTRODUCTION 

Steps on vicinal surfaces exhibit many different insta- 
bilities in the presence of non equilibrium driving forces. 
Of particular interest to us here are the current-induced 
instabilities on Si surfaces that were first discovered by 
Latyshev et alX in 1989. After resistive heating of a vic- 
inal Si(lll) surface with a step-down direct current at 
temperature around 900° C, they observed the formation 
of closely packed step bunches separated by wide step- 
free terraces. The the uniform step train remained stable 
on heating with a step- up current. This instability has 
a mysterious temperature dependence)2i^*^i^ with three 
temperature ranges between 830°C and 1300°C where 
the unstable current direction reverses. 

Furthermore, recent experiments^^ in temperature 
range II (about 1050°C to 1150°C) have shown that after 
heating for several hours with a step-down current, the 
initially uniform steps exhibit a novel wandering insta- 
bility with finite wavelength in-phase sinusoidal undula- 
tions in their positions. When the current is directed at 
an angle to the average step direction, the undulations 
are continuously distorted by the field until finally all the 
peaks point in the direction of the field. 

These instabilities likely arise from a complex inter- 
play between the driven diffusion of adatoms induced 
by the electric field E and their attachment/detachment 
kinetics at steps, which serve as sources and sinks of 
adatoms. (Island formation is not important in the tem- 
perature regimes we consider.) Adatoms are believed to 
acquire a small effective charge z*e, which includes both 
electrostatic and "wind-force" contributions arising from 
scattering of charge carriers, and thus experience a force 
F z*eE that biases their diffusive motioniii Typically 
E has a magnitude of about 5 ~ lOV/cm and z* is of the 



Most theoretical methods are based on a generalization 
of the approach taken by Burton, Cabrera and Frank 
(BCF)fi^ where one considers field driven diffusion of 
adatoms on terraces, with boundary conditions at the 
steps, viewed as line sources and sinks. We will gener- 
ally refer to these generalized BCF models as sharp step 
models. Surface reconstruction typically seen on semi- 
conductor surfaces clearly has important effects on the 
movement of adatoms on terraces and may well affect 
the attachment kinetics at steps. While it is relatively 
simple to take account of reconstruction on terrace dif- 
fusion by changing the diffusion constant, it is much less 
clear how it should be incorporated into the boundary 
conditions at the step edges. 

Many different boundary conditions have been pro- 
posed, incorporating, e.g., asymmetric attachment- 
detachment barriers^^ii^ periphery diffusion along a 
step,^^ permeable steps jiSii^ and field-dependent kinetic 
coefficient)^ and researchers have shown that different 
combinations can give results that can agree with some 
experiments on current-induced step bunching. How- 
ever, a general understanding of the physics leading to 
the sharp step boundary conditions and how they are af- 
fected by reconstruction and the external field is far from 
clear. 

We discuss here a simple model incorporating the key 
physical features of driven diffusion and surface recon- 
struction. It can provide a consistent explanation of 
many experimental results on both Si(lll) and Si(OOl) 
surfaces in terms of a few effective parameters. The 
model also provides a physically suggestive way of inter- 
preting sharp step boundary conditions, showing how the 
effective parameters in continuum models can be related 
to kinetic processes on vicinal surfacesiSi 
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FIG. 1: The upper part of the figure shows a 2D schematic 
view of the vicinal surface composed of different reconstruc- 
tion regions on terraces and near steps, separated by dashed 
fines. In tfiis paper, we assume that the step reconstruction 
with a fixed width s always follows the motion of the step 
(solid line). The lower part of the figure shows a correspond- 
ing ID side view that illustrates our coordinate system. 

In the following, we will set up the basic two-region 
diffusion model and then examine the non-equilibrium 
steady state (NESS) solutions. Step bunching and wan- 
dering regimes are discussed, and combined to provide a 
coherent scenario for the complicated Si(lll) electromi- 
gration experiments. A mapping to the generalized BCF 
model is presented next that sheds light on the sharp step 
boundary conditions. Finally, we study the long time and 
nonlinear behavior of these instabilities and find some in- 
triguing patterns that resemble many features seen in real 
experiments on Si(lll) surfaces. An alternate derivation 
of the basic equations and applications to the different 
instabilities seen in Si(OOl) surfaces is given elsewherei^ 

II. TWO-REGION DIFFUSION MODEL AND 
STEADY STATE SOLUTIONS 

It is well known that the dangling bonds at semicon- 
ductor surfaces quite generally rearrange to form charac- 
teristic surface reconstructions. We expect a different lo- 
cal rearrangement of bonds in the vicinity of a step, which 
itself represents an additional source of dangling bonds. 
Clearly this reconstruction can directly influence surface 
mass transport and hence possible instabilities. Standard 
boundary conditions in the continuum sharp step model 
may include some effects of surface reconstruction in spe- 
cial cases. For example, Liu and Weeks^i interpreted 
electromigration experiments in the lowest temperature 
regime of Si(lll) using attachment/detachment limited 
kinetics, and argued that the attachment barriers could 
arise from a local reconstruction of the dangling bonds 
at a step edge. However, it is not clear how this picture 
should be modified at higher temperatures. 

Steps differ fundamentally from terraces by serving as 



sources and sinks for adatoms. In the classical BCF pic- 
ture it was assumed that the local equilibrium concen- 
tration of adatoms at a step is maintained even in the 
presence of nonequilibrium driving forces. In addition 
the rates of various mass transport processes near steps 
can differ from kinetic processes on terraces, e.g., because 
of differences in local surface reconstructions. The kinetic 
coefficients in generalized BCF models try to take both 
features of steps into account in an effective way. 

Our approach here is to consider a more detailed de- 
scription where both features are treated separately in 
the simplest possible way. We then obtain the relevant 
sharp step boundary condition by an appropriate coarse- 
graining. To that end, we assume that an atomic step 
has sufficient kink sites to maintain a local equilibrium 
concentration of adatoms as in the classical BCF pic- 
ture. Reconstruction is taken into account by assuming 
that the atomic step is surrounded by a step region where 
adatoms undergo effective diffusion with a diffusion con- 
stant Ds that can differ from Dt, the value found on 
terraces. 

Here we use the simplest realization of this idea, where 
the reconstruction is assumed to occur fast relative to 
step motion, so that the step region moves with the 
atomic step and has a fixed width s of a few lattice spac- 
ings at a given temperature. Thus a uniform vicinal sur- 
face can be viewed as an array of repetitive two-region 
units, made up of the nth step region and its neighbor- 
ing lower terrace region. We assume that straight steps 
extend along the y direction and that the step index in- 
creases in the step-down direction, defined as the positive 
X direction, asde schematically shown in Fig. ^ 

The adatoms undergo driven diffusion from the electric 
field. The biased diffusion flux of adatoms with density 
c takes the form: 

F 

Ja = -Da^Ca + Da-, —Ca, (1) 

Kb -I 

where a = {t, s) indicates the terrace or step region and 
Da is the diffusion constant in the corresponding re- 
gion, which here is taken to be isotropic for simplicity. 
We also assume that the effective charge is the same in 
both regions and ignore the small effects of step motion 
on the steady state adatom density field, since the di- 
rect field-induced adatom drift velocity is generally very 
much larger than the net velocity of the steps (driven 
by free sublimation in real experiments) even at high 
temperatures 1^ 

In many studies of step dynamics, because the separa- 
tion of their respective time scales, it suffices to solve the 
diffusion problem with fixed step positions and then bal- 
ance the fluxes locally at a step to determine its motion. 
This is often called the quasi- stationary approximation, 
and it will be adopted throughout this paper. Thus the 
static diffusion problem is simply 

V • J„ = (2) 

in each region, along with continuity of c and J at flxed 
boundaries between terrace and step regions. The normal 
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FIG. 2: Plot of concentration profiles according to Eq.Q with 
model parameters. R — 10 for (i) and (ii), 7? = 0.1 for (iii) 
and (iv); \fa\ = 0.01 in all cases. 



width in the steady state. C is a constant to be deter- 
mined shortly. 

Evidently, it is the interplay between the external elec- 
tric field and changes in the local diffusion rates, charac- 
terized by various combinations of the two parameters / 
and R, that causes the intriguing instabilities. With the 
electric field perpendicular to the step region, altogether 
there are four types of steady state adatom concentration 
profiles with different combinations of parameters / and 
R, as shown in Fig. |21 In the absence of sublimation, 
the concentration profiles we obtain here are completely 
driven by the external field. By taking the limit / — > 
in Eq. @ , one should recover the equilibrium concentra- 
tion (denoted as Ceq) on the entire surface. This fixes the 
constant in Eq. ^ as 



C = c,, {k + s) I {k + Rs) . 



(6) 



Moreover, the constant flux at NESS can be written as 



velocity of the step region is given by mass conservation 
locally at an infinitesimal portion of the step region 



dr [3 s 



(3) 



Here denote diffusion fluxes in the front and back 
terraces respectively and Ac is the difference of the areal 
density of the two phases — the solid phase and the 2D 
adatom gas phase. For simplicity, we take a simple cubic 
lattice, so that Ac « l/fi = a~^, where a is the lattice 
parameter. The last term in Eq. J^J represents the con- 
tribution from diffusion flux in the step region parallel to 
the step, where r denotes the arc length. 

Eqs. define the two-region diffusion model. We 

first consider the NESS solution corresponding to a ID 
uniform step train. In this case, the step normal direc- 
tion coincides with the x direction on terraces, and thus 
parallel or tangential diffusion in the step region plays no 
role here. The NESS concentration profile (denoted by a 
superscript '0') in a two-region unit is easily obtained by 
solving Eq. 10) in both regions subject to continuity of 
concentration and fluxes at the boundaries and is given 
by 



R 



1 - 



(1 - i?) (1 - e 



(4) 



Here 



R = Dt/Ds 



(5) 



is one of the key dimensionless parameters that describes 
the relative diffusion rates of adatoms on terraces and 
in the normal direction of step regions, / = F • x/ksT 
has a dimension of inverse length and characterizes the 
strength of the external field, and It is average terrace 



Jo(/) = DtCeqf 



l + {R-l)s' 



where 



l = lt + s 



(7) 



(8) 



is the distance between the centers of two adjacent step 
regions in a uniform step train. Note that the NESS 
concentration profile of adatoms given by Eq. reduces 
to a constant on the entire surface in presence of the field 
if the diffusion in the normal step direction is the same 
as terrace diffusion, i.e., when R = 1. 



III. STEP BUNCHING AND WANDERING 
INSTABILITIES 

In this section, we study the stability of the NESS so- 
lutions. In particular, the physical origins of both step 
bunching and wandering instabilities are qualitatively 
discussed. 



A. Step Bunching Instability 

A common feature of all NESS profiles shown in Fig. [21 
is that adatom concentration gradients build up in both 
terrace and step regions. Under experimentally relevant 
conditions the field is sufficiently weak that fs < fit ^1 
and linear concentration (or chemical potential) gradi- 
ents form. It is then easy to see that the local equilibrium 
boundary condition c = Ceq in the center of the step re- 
gion holds automatically by symmetry. In the qualitative 
picture of step bunching discussed by Liu and Weeks^ a 
positive terrace concentration gradient (induced in their 
model by a step-down current with an attachment bar- 
rier at a sharp step edge) leads to step bunching. The 
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steady state profile they analyzed leading to step bunch- 
ing in temperature regime I is very similar to case (i) in 
Fig. 121 This corresponds in the two-region model to a 
step-down field with slower diffusion in the step region, 
in agreement with an intuitive picture of a step barrier. 

Moreover, it is clear that profile (iv) is qualitatively the 
same as (i). Hence we expect that the steady state (iv), 
corresponding to faster diffusion in the step region with 
a step-up field, also undergoes a bunching instability. A 
hopping model with these features was studied by Suga et 
alm^ by computer simulations, and indeed they observed 
a bunching instability. 

To understand the bunching of straight steps it is use- 
ful to consider a ID version of Eq. © : 



Vn = Vl [Jo {In-l) - Jo (L)] , 



(9) 



where the ID flux Jg as given by Eq. Q now depends 
on the local terrace widths. Consider a small devia- 
tion Sxn = e„e'^'^* for nth step from the NESS, where 
En = ee™*^ and </> is the phase between neighboring steps. 
Then the step will move as a result of the unbalanced 
fluxes induced by changing width of the terrace in front 
= l + £^ (e*</> - 1) and back l^^i = I + {l - e'"^). 
The amplification rate uji is given by ui — and 
substituting into Eq. © gives 



UJi 



C0S( 



dl 

90n n 



(10) 



(i?- l)s]2 



COS( 



Clearly, step bunching occurs when / (i? — 1) > 0, cor- 
responding to two different regimes discussed above, and 
in both cases the most unstable mode is a step pairing 
instability with = tt. 



B. Step Wandering Instability 

The ID NESS concentration profiles also provide im- 
portant insights into step wandering, which is essentially 
a 2D phenomenon. It is clear that the concentration gra- 
dient on the terraces in cases (i) and (iv) can drive a step 
wandering instability. The monotonically increasing ter- 
race chemical potential tends to make a forward bulging 
part of a step move even faster, as was first demonstrated 
for vicinal surfaces by Bales and Zangwill.15- This is 
the essence of the classic MuUins-Sekerka instabilityi2Si2i 
However, as shown above, these same profiles lead to ID 
step bunching, which tends to suppress the wandering in- 
stability. Moreover, this mechanism cannot explain the 
behavior in regime H of Si(lll) where wandering and 
bunching occur for different current directions. 

The fact that this step wandering cannot be of the 
Mullins-Sekerka type driven by terrace gradients suggests 
that it may originate from mass transport in the step re- 
gion. Let us focus on a single 2D step region, as in Fig. 13 
In this case, it is convenient to describe the step region 




FIG. 3: A geometrical view of a single wandering step region. 
The dashed arrows inside the step region schematically shows 
the driven flux that is parallel to the step for a step-down {x 
direction) field. The lower right corner shows the case when 
the field is at an angle ip off the a::-axis. 



using curvilinear coordinates set up by the local normal 
and tangential directions of the step. For a long wave- 
length step fluctuation with wavenumber q there exists 
a nonzero component of the field in the tangential direc- 
tion, which induces a driven flux along the step propor- 
tional to fq^. For a step-down field (/ > 0), this driven 
flux is destabilizing since it tends to transport mass from 
"valleys" to forward-bulging "hills" . On the other hand, 
the stabilizing flux due to the curvature relaxation is pro- 
portional to Fg^, where F is an effective capillary length 
in the step region. The competition between these two 
terms results in a finite wavelength linear instability, oc- 
curring on a length scale of order ^, where 



(11) 



In principle this new wandering instability could arise 
in cases (i) and (iii) of Fig. [21 where there is a step-down 
field. However step bunching also occurs for case (i). 
Only case (iii) with / > and faster diffusion in the step 
region (R < 1) is free of step bunching, and thus capable 
of explaining experiments in Regime II of Si(lll). In the 
next section we show that these qualitative conclusions 
are in agreement with a more detailed analysis based on a 
mapping of the two-region model to an equivalent sharp 
step model. 



IV. MAPPING TO A GENERALIZED BCF 
MODEL 

In this section we show how the two-region model can 
be used to generate the appropriate sharp step boundary 
conditions by a mapping to a generalized BCF model. 

The general continuum boundary condition in the 
sharp-step model assumes small deviations from local 
equilibrium and introduces linear kinetic coefficients k± 
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FIG. 4: Shown is a highly exaggerated profile for a downhill 
force and slower diffusion in the step region. Also illustrated 
with the dashed-dot line is the extrapolation of the terrace 
profile to the center of the step region, thus determining the 
parameter cf in Eq. 11211 . The lower part of the figure gives a 
side view of sharp equilibrium steps and their associated step 
regions. 



to relate (or c[), the limiting lower (or upper) ter- 
race adatom density at the step edge, to the associated 
terrace adatom flux into the step. To linear order in the 
field this gives rise to the standard sharp step boundary 
condition: 



±Dt [Vct - fCt]^ ^ k{ct- Ceq) 



(12) 



Here k is the corresponding sharp step kinetic coefficient, 
which is symmetric in this case. 

A natural way of relating the NESS solutions of the 
two-region model to those of sharp step model is to ex- 
trapolate the terrace concentration profile to the center 
of the step region. This corresponds to a physical coarse- 
graining where the step region has negligible width when 
compared to the terrace widths. The use of extrapola- 
tion to relate the parameters in discrete and continuum 
models is well known in other interface applications^ 
We use Eq. 10} to evaluate the gradient, and identify cf 
as the extrapolated value of terrace concentration at the 
atomic step, as illustrated in Fig. 01 Substituting into 
Eq. p2|) . to lowest order in the field we find that 



d = 



Df 



1 



(i?- l)s. 



(13) 



Note that the terrace width I in the sharp step model 
is naturally related to the two-region width It hy I = 
It + s, as in Eq. |(SJ. Here d is often referred to as the 
attachment-detachment length. 

Equation H13|l gives a mapping of the parameters in 
the simplest two-region model to those of a generalized 
BCF model. When R > 1 (faster diffusion in the terrace 
region), k is positive, which leads to a bunching instabil- 
ity for a step-down current. When R = 1 (the diffusion 
rate is the same in both regions), k goes to infinity, which 



forces cf in Eq. (|12(l to equal Ceq , corresponding to local 
equilibrium with no instability. When R < 1 (diffusion 
is faster in step regions than in terrace regions), k be- 
comes negative, which leads to step bunching by a step- 
up current together with step wandering by a step-down 
current. 

The possibility of a negative kinetic coefficient, or 
equivalently a negative d, was first suggested in the work 
of Politi and Villain,— though with no derivation or dis- 
cussion of any physical consequences. Note that even 
though the derivation given here considers a terrace con- 
centration profile obtained by electromigration, Eq. H13|l 
is a general result that is independent of the field. In a re- 
lated workjSS we derive sharp step boundary conditions 
by considering a discrete hopping model with different 
hopping rates in two regions but without the field, and 
again obtain Eq. Ijl^fl . 



V. LINEAR STABILITY ANALYSIS 

With the mapping defined by Eq. l(T^ . the hnear sta- 
bility analysis can be performed using a standard sharp 
step model, with parameters obtained from the physi- 
cally suggestive two-region model. The general result is 
presented in the appendix, and here we concentrate on 
the resulting stability in the weak field (/Z <C 1) and long 
wavelength {ql <^ 1) limit. The real part of the stability 
function can be written as 



UJr UJl if, (t>) + t^2 [q, /, ' 



where 



= flDtcl 



eq 



{I + 2d)' 



(1 



COS ( 



(14) 



(15) 



and 



L02 = ^Dtc^^q^ \ -V 



+f 



2(1 — cos ( 

l + 2d 
2dl s 



R 



l + 2d R 



(16) 

ioi characterizes the ID instability and thus is inde- 
pendent of q. The bunching instability occurs for df > 
with most unstable mode giving step pairing with (j) = n. 
Note that Eq. (HHJl is identical to Eq. when Eq. ifT^ 
is used. 

UJ2 characterizes 2D wandering instabilities with re- 
spect to perturbations of wavenumber q. The first term 
on the right hand side is stabilizing, and has its minimum 
value for ip = 0, where it is proportional to Fq* and all 
the steps wander in phase. 

The second term, proportional to the field, contains 
two destabilizing contributions. The first contribution, 
proportional to Dtdfq^, describes a MuUins-Sekerka or 
Bales-Zangwill instability induced by the terrace concen- 
tration gradient that can occur when df > 0. 
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TABLE I: Linear Stability Results 





d>0 {R>1) 


d < (7? < 1) 




Bunching with maximum mode (j) — n 
Wandering with maximum mode = 


Wandering with maximum mode (p — 


/<o 


Linearly stable 


Bunching with maximum mode (p = n 



The second contribution, proportional to Dgsfq^, rep- 
resents an alternative mechanism for step wandering in- 
duced by field-driven periphery diffusion along the step. 
When d > 0, both mechanisms operate with a step- 
down current, while the step-up case is completely stable. 
When d < 0, the second mechanism can produce wander- 
ing with a step-down current, while bunching occurs for 
a step up current, as was discussed earlier in Sec. (|fff Bll . 
These stability results are summarized in Table 



where i? = f. Conceivably, such a transition could 
happen again at higher temperatures, since only small 
changes in the relative diffusion rates can take the fun- 
damental parameter R from less than to greater than 
unity and vice versa. This scenario provides a consistent 
interpretation of experiments in the second temperature 
regime and suggests more generally why there could be 
such a complicated temperature dependence. 



VI. IMPLICATIONS FOR SI SURFACES 



VII. NONLINEAR EVOLUTION OF 
CURRENT-INDUCED INSTABILITIES 



Thus far, both step bunching and wandering instabil- 
ities have been analyzed in general terms based on the 
simple idea of two-region diffusion. Now we examine the 
implications for vicinal Si(lll) surfaces. If we assume for 
concreteness that reconstruction is generally associated 
with slower adatom diffusion, we can give a qualitatively 
reasonable scenario that can account for many features 
of the electromigration experiments observed on Si(lll). 

In temperature range I, we assume there exists recon- 
struction in both step and terrace regions. Consistent 
with the analysis of Liu and Weeks, we assume that 
at low temperature the adatom diffusion in the recon- 
structed step region is slower than in the terrace region, 
i.e. R> 1, corresponding to cases (i) and (ii) in Fig. [3 A 
step-down current induces both step bunching and step 
wandering of MuUins-Sekerka type. However the wan- 
dering is likely suppressed by the bunching instability. A 
step-up current produces a stable uniform step train. 

At higher temperatures, we expect reconstruction in 
step region could have a more fragile structure when 
compared to that in the terrace region since step atoms 
have more dangling bonds. Thus there could exist an in- 
termediate temperature range where because of changes 
in the step reconstruction, diffusion is faster in the step 
region than on terraces, i.e. R < I, corresponding to 
cases (iii) and (iv) in Fig. |21 The uniform step train 
now exhibits bunching with a step- up current. Wander- 
ing occurs with a step-down current, induced by driven 
diffusion parallel to the step. In particular, if we sub- 
stitute in Eq. Ijlll) the latest experimental values for 
the step stiffness)^ /3 — 16.3meV/ A, and for the ef- 
fective charge)24 z* = 0.13, and use a typical electric 
field strength oi E = 7V/cm, the resulting wavelength 
is roughly given by A ~ 27r^ ~ 5/im, comparable with 
experimental values^^ of 6 — dfim. 

In this picture, the transition between different tem- 
perature regimes is associated with local equilibrium 



A. Velocity Function Formalism 

To calculate the long time morphology of vicinal sur- 
faces, effective equations relating the velocity of a step to 
the local terrace widths have proved to be very usefuli^ 
A simple example of such a velocity function is given by 
Eq. Q. The extended velocity function formalisn&^ 
takes into account also the capillarity of steps (line ten- 
sion effects) as well as step repulsions, which are needed 
to prevent step overhangs as the initial instabilities grow. 
Here we also incorporate a periphery diffusion term, the 
sharp step analogue of the parallel diffusion flux in the 
two-region model. Thus the general form of the velocity 
function can be written as: 

«n(2/) = /+ (^n;Mn:Mn+l) + /- (^n- 1 i 1 , Mn) - dr^s 

(17) 

where Z„ (y) is the local width of terrace n that is in front 
of step n and fin{y) is the local chemical potential of the 
step n. 

The velocity functions f± contains contributions both 
from driven fluxes on the two neighboring terraces given 
by the sharp step equivalence of Eq. (|5J), and equilibrium 
relaxation terms that can be calculated in terms of the 
step edge chemical potentials The /i„ take account 

of both capillary effects for an individual step (using a 
linear approximation for the curvature) and the effects 
of nearest neighbor step interactions as described earlier. 
See Refs. (321 and 13^ for detailed expressions for f± and 

Numerically integrating Eq. I|17|) . we find step bunch- 
ing patterns for two parameter regimes (i) / > 0, i? > 1 
and (ii) / < 0, i? < 1, in agreement with predictions of 
linear stability analysis. The bunching patterns in these 
two regimes are qualitatively similar, as shown in Fig. \S[ 
In both cases, step bunches form and grow. In between 
the step bunches there are crossing steps traveling from 
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FIG. 5: A uniform step train composed of 30 steps with spac- 
ing of i = 10 forms step bunches at later times both for (i) 
/ > 0, i? > 1 and (ii) / < 0, 7? < 1. 




FIG. 6: A uniform step train comprised of 5 steps with spacing 
I — 5 forms in-phase wandering patterns at later times for 
f > 0, R < 1. Notice there are some defects in the pattern 
because the wandering wavelength is incommensurate with 
the finite size of our system in the y-direction. 



one bunch to the other. 

In-phase step wandering is also given by Eq. IjlTI) in the 
regime / > 0, i? < 1, as suggested by the previous lin- 
ear stability analysis. Typical wandering patterns with 
model parameters are shown in Fig. |3 Even though this 
is known to be a linear instability, numerically we observe 
that it acts very like a nucleation process. The steps fluc- 
tuate randomly as if the surface were completely stable 
until a sinusoidal perturbation of the right wavelength 
forms. Once formed, these small scale sinusoidal waves 
propagate through effective "pulling" by capillary effects 
in the lateral direction and by step repulsions in the nor- 
mal direction, until the entire surface is covered. This 
is qualitatively consistent with experimental findings on 
Si(lll)^ 



B. Evolution of Step Wandering in a Geometric 
Representation 

Although Eq. (|17|l has captured many physical fea- 
tures, it uses a linearized curvature approximation and 
cannot be trusted when the step curvature becomes large. 
Recent experiments show a continuous distortion of the 
sinusoidal wandering wave by a field directed at an angle 
to the step normal. We treat this problem here using a 
geometrical representation^^*^ of the step, where a sin- 
gle curve is parameterized by intrinsic properties like its 
arc length r and curvature n. 

It suffices to concentrate on a single step, since step 
wandering occurs in phase. Consider a geometric repre- 
sentation of our step region with constant width s, as in 
Fig. O The morphology of the step region is specified 
by the position vector x {t, t) of the atomic step in the 
middle, where r can represents the arc length measured 
from an arbitrary origin. To follow x {t, r) at a later time 
we need to know the velocity of the curve 



9x 
'dt 



Vnfl + Vrf, 



(18) 



where n and f denote normal and tangential directions 
as before. 

A general treatment of time-dependent curvilinear 
coordinates^*', shows the equation of motion for the curve 
is 



Ok 
'dt 



Ok 

OT 



(19) 



which is subjected to the nonlocal metric constraint 
dr 



dt 



Vr (r) — Ut- (r = 0) + / VnKds' 



(20) 



Interpreting r as the arc length is arbitrary and other pa- 
rameterizations can be used, since only the normal veloc- 
ity of the curve is physically relevant. Following previous 
workers,'^- we take advantage of this "gauge freedom" and 
choose the orthogonal gauge, where r is chosen at each 
instant of time so that the interface velocity has only a 
normal component (vt- — 0). 

Now, we need to determine the normal velocity along 
the step. For simplicity, we will neglect contributions 
from the terrace diffusion field as well as from the nor- 
mal diffusion field in the step region, since it has already 
been shown that the wandering instability we are inter- 
ested in is induced by the biased diffusion parallel to the 
step. In the quasi-stationary limit, the diffusion field 
inside the step region is stationary for any given step 
position. To a good approximation, it can be taken as 
Cs — (1 -t- Tk), where = c^^s is the adatom density 
per unit step length for straight steps. 

Next we consider the time rate of change of the 
adatoms contained in an element of the step region with 
an infinitesimal length 5t that moves with velocity w„ 
as in Fig. 13 This balance contains contributions from 
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FIG. 7: Step evolution under a perpendicular electric field 
(a) At t = 160, a linear instability develops; (b) At t = 170, 
asymmetry between the peaks and valleys creates a periodic 
cellular structure; (c) At t = 190, the cellular shape is pre- 
served but it grows in amplitude. 



FIG. 8: Step evolution when the electric field is at an angle 
if = 7r/4 from the x-axis: (a) t — 300, the initial instability 
induced by the normal component of the field; (b) t — 315, 
the peaks have begun to turn; (c) t = 330, all the peaks align 
with the direction of the field. 



the motion of the step, and from the divergence of the 
flux parallel to the step. The latter accounts for diffusion 
driven both by the field and by chemical potential vari- 
ations arising from changes in step curvature. We thus 
have 



1 ^'^^^^ 



ft ^VnSr — Dgdr [fcs sin {(p — 6)] 5t 

2„ x_ (21) 



Using the exact geometrical relation [d{5T) /dt]^ = 
VnKSr, which can be understood physically as the rate 
at which the arc length St on a circle of radius 
changes if the circle grows only radially at rate u„, Eq. 
(|21|1 reduces to the following form 

Vn [H-f7c°(l + rK) k] 



/cos {ip~0)il + TK)K 

-f sin (<y9 - e) TdrK + Td^K. 

(22) 

Combining Eq. if^ with Eqs. ifT^ and (HOI yields a com- 
plete description of the dynamics of a single step region 
in the presence of an electric field at an angle ip off the 
X-axis. 

We first consider the special case (p = where the ex- 
ternal field is perpendicular to the average step direction 
(the y -axis). In Fig.[7| we show three step configurations 
evolving from a straight step with a small perturbation 
in the middle. The linear wandering instability develops 
first as shown in Fig.[7fa), then gradually changes into a 
cellular shape with the wavelength selected by the linear 
instabihty, as illustrated in Fig. EKb). At later stages, 
the cellular shape grows without significant distortion or 
overlap, as shown in Fig. [Tfc). Notice that indeed we 
observe numerically a long time period before the linear 
instability is significant. 




40 60 

(b) 
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FIG. 9: Comparison of the step evolution as the angle ip in- 
creases: (a) t — 230, (p = 7r/6; (b) t = 330, — tt/4; (c) 
t = 640, p = 7r/3. 



In Fig. |S1 we show configurations of the system with 
ip = 7r/4. Fig. |S1 suggests that the linear instability is in- 
duced by the perpendicular component of the field. How- 
ever, as the magnitude of the instability grows, the peaks 
turn gradually until they are aligned with the direction 
of the field. We see the same peak turning process when 
the angle ip is varied while keeping / constant. How- 
ever, since the perpendicular component decreases with 
increasing ip, both the wavelength selected by the ini- 
tial instability as in Eq. Ullfl and the time period before 
it forms increases monotonically with ip. The numerical 
results for three particular angles are shown in Fig. O 

To provide a more qualitative understanding of the 
pattern formation process, we neglect the higher order 
terms in k in Eq. H22|) . To linear order in k, Eq. H22|) 
becomes 

fK cos {ip-9)-f sin {ip ~ &) drK + Vd^K. (23) 
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FIG. 10: A study of the asymmetry of the cellular patterns; 
(a) t — 180, a snapshot of the system given by Eq. 
Note the close agreement with Fig. |7[c). This shows that 
the simphfed Eq. I|28|l with only terms linear in k captures 
most features of Eq. 1221 : (b) t — 180, a snapshot of a model 
equation where the term ~ is left out of Eq. I{'24^ . Clearly 
this term is mainly responsible for the asymmetric shape in 
(a). 



In particular, for ip = 



no. 



fKC0s9 + / sin 9drK + d^K. 



(24) 



In the usual MuUins-Sekerka instability k alone ap- 
pears in the first term. Here however we have kcos0, re- 
sulting from field driven diffusion inside the step region. 
The extra cos 9 term brings in a field induced anisotropy 
that makes the peaks and valleys of a perturbation prefer- 
ably grow rather than the sides. This stabilizes cellular 
structures. This anisotropy will keep the tip unsplit, and 
it provides a cut off as the sides become nearly vertical. 
Thus the cellular shapes formed under the influence of 
the external field do not emit side branches, in contrast 
to most systems that undergo a MuUins-Sekerka instabil- 
ity. 

The second term in Eq. H23() is a flux induced by — k 
that effectively transports mass from the bottom to the 
top of a bulge and is responsible for the asymmetric shape 
of the peaks and valleys, as is illustrated in Fig. ^| 

Although Eq. (|23|l is linear in the curvature, k itself 
is a highly nonlinear function of the deviation from a 
straight step. The early evolution is governed by the 
following linearized equation 



1 



dx d'^x J. . d^x „u'x 

^- = -/cos^^-/sm^^-r^. (25) 

The above equation is unstable when f cos if > 0, sug- 
gesting that the wavelength selection is determined by 
the perpendicular component of the field. For — Q, 
perturbations with wavenumber ~ l/(-\/20 are max- 
imally amplified. For < (/s < 7r/2, the most unstable 
wavenumber selected by the linear instability is decreased 
by a factor of ^cos Lp, i.e., q^p = go-y/cos Lp. 



As the instability grows, the field induced anisotropy 
characterized by the factor cos,{if — 9) becomes more sig- 
nificant. As in the Lp = Q case above, the anisotropy 
makes the initial sinusoidal wave grows preferably in the 
direction where cos{(f — 9) in Eq. H23I) attains its mini- 
mum. Thus the wave will be continuously distorted until 
the peaks point toward the field direction, and subse- 
quently only the magnitude of the pattern grows. 



VIII. CONCLUSION 

In this paper we have studied a physically suggestive 
two-region diffusion model. The basic idea is to consider 
different hopping rates associated with different recon- 
struction and bonding in the terrace and step regions. 
The resulting steady state profiles provide important in- 
sight into the physical origins of both step bunching and 
wandering instabilities. Step bunching is induced by pos- 
itive chemical potential gradients on terraces that are 
essentially determined by the sign of /(i? — 1). We ar- 
gue that step wandering in Si(lll) does not arise from 
the well known MuUins-Sekerka instability. Rather, it 
is induced by driven diffusion along the step edge under 
the influence of a step-down force, and only becomes sig- 
nificant when step bunching is absent, which requires a 
negative kinetic coefficient. 

We also carried out a mapping from the two-region 
model to a sharp step model using a simple extrapolation 
procedure. The result connects the kinetic coefficients in 
sharp step models to relative diffusion rates in terrace 
and step regions. In particular, the lowest order result 
shows that the kinetic coefficients are independent of the 
driving field, in contrast to earlier suggestionsiSS 

A coherent scenario for Si(lll) electromigration is pro- 
posed based on the stability analysis of the model. In 
particular, the mysterious second temperature regime is 
interpreted using a negative kinetic coefficient. This al- 
lows the step wandering that generally occurs with a step- 
down force to be separated from step bunching. The 
transition between different temperature regimes is gov- 
erned by the relative diffusitivity in the terrace and step 
regions. Other theories can predict a reversal of step 
bunching arising from a change in step transparencjii£ii^ 
or from a change of sign of the effective chargei^ How- 
ever, neither approach can give a consistent treatment 
for step wandering. 

The long time evolution of the step instabilities was 
calculated by numerical integration of a set of equations 
based on the standard velocity function formalism with 
the addition of a periphery diffusion term. The linear 
instabilities are recovered at short times and interesting 
2D pattern formation is see at longer times in qualitative 
agreement with experiment. 

We also showed that a geometric representation of the 
step provides a simple way to describe the nonlinear evo- 
lution of step wandering patterns with large curvatures. 
The resulting cellular patterns when the driving field is at 
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an angle to the step shows significant step "overhangs", 
which can not be captured by standard multi-scale ex- 
pansion methodsi^2iili 

The two-region model can also be modified to explain 
many features of the very different step bunching behav- 
ior seen on Si(OOl) surfacesi^ Thus it provides a sim- 
ple and unified perspective that can shed light on both 
general properties of current-induced step bunching and 
wandering instabilities and their specific manifestations 
on Si surfaces. 



Consider a 2D perturbation on the step profile in the 
form Sxn (y, t) = Xn {y, t) — = e^e'^*"'""''' -I- c.c, where 
is the step position for ID steady state and £„ is the 
ID perturbation previously defined. In general lo can be 
complex, i.e. to = + i^i, but we are only interested 
in the real part Wr whose sign determines the instability. 
The calculation follows standard methods, and the result 
can cast in the familiar Bales and Zangwill's form^'"': 

= -Tq^h (q, /, <f>) + fg (q, /, 0) . (A.5) 
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APPENDIX: LINEAR STABILITY ANALYSIS IN 
A SHARP INTERFACE MODEL 

A complete 2D stability analysis in a generalized BCF 
model is performed in this section, with boundary con- 
ditions dictated by mapping from the two-region model. 
Using the quasi-stationary approximation, we first solve 
for the static concentration field ct on the terrace as given 
by 



AF 



^ct = 0, 



(A.l) 



subject to the general linear kinetics boundary condition 
at the sharp step: 

±Dt [Vet - /ct]± ■h^k[ct- (1 + Tk)] ^ . (A.2) 

Here n is the curvature, defined to be positive for a circle. 
The normal step velocity is determined by balancing the 
fiuxes locally at the step 



(A.3) 



Here Js is the periphery fiux of the mobile atoms along 
the interface, which can be viewed as the coarse-grained 
contribution from the parallel diffusion in the two-region 
model. In general, Js takes the form 



F • f 



(A.4) 



where Cg — c^^s gives the effective number of ledge atoms 
per unit step length. 



Here the stabilizing piece h (q, /, (f>) is given by 

h{qj,(j)) _ 2A[cosh(A/)-cos(/>cosh(/V2)] + 2dq2sinh(A/) 

nDtc° 



eg 



V 



a o 



(A.6) 



and the destabilizing piece g (g, /, 4>) is 
9{qJA) 2df 



[df (e/' + 1) + e/' -1]V 
X {2A [cosh(A/)-e/'/2cos(^] 
-|-2dg2e/'/2 cosh if 1/2) sinh {XI) 
+ sinh(/V2) (A+e^-' - A_e^+')} 



where 



V = 2dA cosh (AO + (l + d'^q^) sinh {XI) , 



(A.7) 
(A.8) 



A = + V/2 and A± = //2 ± A. 

It is easy to see that h {q, /, (j>) is positive definite; thus 
the first term in Ea. (|A.5|l is always stabilizing. In partic- 
ular, we obtain the results for the equilibrium relaxation 
by taking the limit / — ^ 



^0 {q, 



2/2? [cosh {ql) — cos 4> + dq sinh {ql)] 
^ \2dgcosh((j0 + (l-|-d2(j2)sinh {ql) 



s 2\ 



(A.9) 



The two terms in the curly brackets account for relax- 
ation through terrace diffusion and periphery diffusion 
respectively. The second term in Eq. IjA.Sp is completely 
driven by the field, and it vanishes identically as / 
since limj„>o 9 {q, f, 0) is finite. 

For current-induced instabilities that are of interest of 
this paper, we can take the weak field {fl ^ 1) and long 
wavelength {ql <^ 1) limit. In this limit, the stability 
functions to linear order in the field are given in Eqs. 
(|TH-(|TC|). 
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